function [ P ] = getPfromPhy( Phy, U_vec, Us_vec, nu, h, NIt, NJt, NIe, NJe, NIut, NJut, NIvt, NJvt, NIue, NJue, NIve, NJve, dimU, dimV, D )
%GETPFROMPHY
%

diffU_vec = U_vec - Us_vec;
DdiffU_vec = D * diffU_vec;

dimP = numel(Phy);

RHS = zeros(dimP,1);
M = zeros(dimP,dimP);

for i = 1 : NIt;
    for j = 1 : NJt;
        if (i<=NIe && j>NJe)
            continue;
        end;
        
        % inedxes calculation
        indexPhy_cen = PmapOpt(i,j, NIe, NJe, NIt);
        
        % Left boundary
        if (i == 1)
            indexPhy_left  = indexPhy_cen;
        else
            if (i==NIe+1 && j>NJe)
                indexPhy_left  = indexPhy_cen;
            else
                indexPhy_left = PmapOpt(i-1,j, NIe, NJe, NIt);
            end;
        end;
        % Right boundary
        if (i==NIt)
            indexPhy_right = indexPhy_cen;
        else
            indexPhy_right = PmapOpt(i+1,j, NIe, NJe, NIt);
        end;
        
        % bottom boundary
        if (j==1)
            indexPhy_bot = indexPhy_cen;
        else
            indexPhy_bot = PmapOpt(i,j-1, NIe, NJe, NIt);
        end;
        
        % top boundary
        if (j==NJt)
            indexPhy_top = indexPhy_cen;
        else
            if (j==NJe && i<=NIe)
                indexPhy_top = indexPhy_cen;
            else
                indexPhy_top = PmapOpt(i,j+1, NIe, NJe, NIt);
            end;            
        end;

        
        
        %Dx = getDx(i, j, diffU_vec(1:dimU), diffU_vec(dimU+1:end), NJue, NIue, NJut, h);
        %Dy = getDx(i, j, diffU_vec(1:dimU), diffU_vec(dimU+1:end), NJue, NIue, NJut, h);
        indexU_left = Umap (i,j);
        indexU_right = Umap (i+1,j);
        indexV_bot = dimU+Vmap (i,j);
        indexV_top = dimU+Vmap (i,j+1);
        
        Dx = (DdiffU_vec(indexU_left) + DdiffU_vec(indexU_right) ) / 2;
        Dy = (DdiffU_vec(indexV_top) + DdiffU_vec(indexV_bot) ) / 2;
        
        % mount Right hand side vecotr.
        RHS(indexPhy_cen) = Phy(indexPhy_left) - Phy(indexPhy_right) + Phy(indexPhy_top) - Phy(indexPhy_bot) +...
            nu*h*( Dx + Dy);
        
        M(indexPhy_cen, indexPhy_left) =  M(indexPhy_cen, indexPhy_left)-1;
        M(indexPhy_cen, indexPhy_right) = M(indexPhy_cen, indexPhy_right)+1;
        M(indexPhy_cen, indexPhy_top) =  M(indexPhy_cen, indexPhy_top)+ 1;
        M(indexPhy_cen, indexPhy_bot) = M(indexPhy_cen, indexPhy_bot)-1;
        
    end;
end;

% M(PmapOpt(NIt,1, NIe, NJe, NIt), :) = 0;
% M(PmapOpt(NIt,1, NIe, NJe, NIt), PmapOpt(NIt,1, NIe, NJe, NIt)) = 1;
% RHS(PmapOpt(NIt,1, NIe, NJe, NIt)) = 0;

P = M \ RHS;

end

